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AN IMPLICIT RUNGE-KUTTA METHOD FOR GENERAL SECOND ORDER ORDINARY 
DIFFERENTIAL EQUATIONS 


SA AGAM, YA YAHAYA AND SC OSUALA 


Abstract. This paper focuses on the derivation of a fully implicit Sixth order Runge-kutta type method with 
error estimation formula for the solution of general second order ordinary differential equations (ODEs). We 
define a transformation on the set of coefficients (Butcher coefficients Tableau) for first order differential 
equations to generate a new coefficient tableau for direct solution of second order ODEs. The scheme is 
simple, A-stable, highly efficient and has low implementation cost. Some problems are used as experimental 
examples. 


1.0 Introduction 

Development of Runge-kutta method for direct solution of second order Odes started with 
Nystrom [1]. He extended the fourth order explicit Runge-kutta method for first order 
ODEs to second order differential equations. Earlier researchers include Huta [2], Felhberg 
[3], Chawla et al [4], Sharp et al [5], Imoni et al [6] Coper et al [7], Filippi et al[8], Hairer [9] 
and Agam et al [10] etc. All these methods mentioned above are explicit type. Explicit 
Runge-kutta methods are unsuitable for solution of Stiff problems because their region of 
absolute stability is small, in particular they are bounded. 

The instability of explicit methods motivates the development of implicit methods. Implicit 
Runge-kutta methods were earlier proposed by Kuntzmann [11] and Butcher [12]. Their 
proposed methods are based on Gauss quadrature. The remarkable thing about these 
methods are that their order p — 2s, for an S-stage scheme and are all A-stable [13] 

The Gauss-quadrature method is a generalization of integrals in the interval (-1,1). The 
Gauss-quadrature methods mentioned above are for first order differential equations. Thus 
there is need to develop parallel Runge-kutta methods for second and higher orders ODEs. 
In this paper we have tried to address these problems by deriving Runge-kutta method 
based on Gauss-quadrature for second order ODEs. 


11  Preliminaries/ Basic definition 

1.2 Definition: Let (G4 *) and (G, .) be two set G4 and G, with binary operations (*) and 
(-) respectively. A transformation T:G, > G, is a linear monomorphism if for any two 
points a4, a2 € G4 
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T(a, * a2) = T(a4,) . T(a;) and T(a,) = T(a2) > a, = a; 

Remark: A monomorphism operators preserve the algebraic structure of the domain into 
its co domain. 

1.2 A Butcher Tableau of coefficients based on Gauss Lengendre quadrature [14] of 
order 6, for solving first order ordinary differential equation is given below 

Table1 Butcher Coefficients table of order 6 for first order differential equations 











C A 
1 15 5 2 15 5 i15 
6-3) s (5-55) (35-50) 
1 5 15 2 5 15 
Gs ui) Gui) 
1 i15 5 15 2 VI5 5 
GeF) P ($435. 36 
bT 5 4 5 
18 9 18 
C A 











bt (b,, bs, bs) 


where A = (aij) 

C = (€4, C2, .... C5) 

bT = (by, bo, ...., bs) (1.0) 
1.4 Error of a Runge-kutta method is defined as 

y(x;) — y; = (Exact error) 

Where y(xj) is the exact solution at x = x; and y; is the RK solution at x = xj. 


2.0 Derivation method 
We consider the general second order ODEs 


azy” + a4y' + aoy = g(x, y), Y' (xo) = y'o (2.01) 
where a;i = 0,1,2 are scalar constants or functions. We then rewrite 2.01 in the form 
y" = fay), yo) = yo y 9 = y', (2.02) 


= f@),v=yy") 
We use the coefficient tableau (1.0) € Ras our basis and domain of transformation T 
defined as follows. 
3 3 
vi=|x+cih,y + > ajT(v;),y' + 2: a;jT (vj) | e R? 


j=1 j=1 
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3 3 3 


T(w)2T|x-ccjh,y-^ 2: ajT(v;), y' +) aiT (v) |-h|y'* 2. aijT (vj) 


j=l j=l j=1 
and 
3 3 
T'(vj) 2 hf| x+cih,y + y ajT(v;),y' + 2 a; T (vj) | - h mj 
j=1 j=1 
3 3 
ie m =f|x+ch,y+ 2. a;T(v;), y' +) aijT (vj) 
j=1 j=1 
Theorem 1 


| 


The Transformation T: R? — R in (2.03) is a well defined monomorphism. 


Proof: 


let u, v € R defined by 
3 3 


j=1 j=1 
3 3 
j=1 j=1 


by the definition of T on R? ,we have 


3 3 
TU+V)=h| y, 2: a;T (uj) +y', + > ajT'(v;) 
j=1 j=1 
3 3 
=h y'i +) ayT (u) +h 35 +) agT (vj) 
j=1 j=1 
= T(U)+T(V) 


Hence T is a homomorphism. 
Now we show that T is one to one 
Letu, v € R? with T(u) = T(v), then by definition of T, we have 
3 3 
h y» + 2 a;T (uj) —h V2 + > aj;T (vj) 
j=1 j=1 
Since T(u) = T(v) then T'(u;) = T'(vj) 


Hence y, — y; and so 


3 
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(2.03) 


(2.04) 


3 3 3 
x+ cih, yi + > aT (u), y +) sr) = |- + cih, V2 + > a;T(v;), y'z +) aron) 


j=1 j=1 j=1 


j=1 
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Hence U = V 

Thus T is one to one > monomorphism 

Remark: By implication the domain of T and the image of T have the same algebraic 
structure. 

We now use this transformation to generate new coefficients and Butcher tableau for 
direct solution of (2.02). 


Using this transformation we have 
3 3 
T(r) - h| y+) ayT(v)) | h| y * 9 a, hm; 
j=1 j=1 
= h(y' + ahm; + a44 hm; + a,3hm3) , where (aij € Table 1) 


5 2 v15 5 V15 


36 9 15 
Similarly, 
3 3 
T(v;) =T x+ czh,y+ ) azjT(v j) y! +) a5;T (vj) 
j=1 j=1 


= h( y' + aj íhim, a; hm, + az3hm3) 


sp TE L ee r Saruy 
Jr aet o umor h ae yg 


T(v3) = h(y' + a34T' (vi) + aa?T' (v5) + assT' (v3)) 
= (y + (+ 2 +e) hom e ($e XE) mes : shims) 


36 30 9 15 36 
By the definition of T, from (2.03) 
3 3 
m, = T'(V;) =f|x+c;h,y + > a;T(v;), y' + > aijT (vj) 
j=1 j=1 


Now substituting for T(V;) from (2.03 ), we have 
m, = (fx + cl, y + aT (v1) + a4QT(v;) + a3T(v3), y' + Ay T' (v1) + Ay2T'(V2) 


+ ai3T'(v3)) 
3 


1 v15 5 . 2 V15 
m =f (s«(- X hv * x^ y * 9 ayT(v;) +(5- =A y * 9 ayT(v;) 


j=l j=l 
Fe seme ED PEE: T puce MEM MEE 
a6 30. y+) ajT(vj) Y ach t (5 7 Ge ma tas 7739 J hms J. 


j=1 
(T'(vi) = hm;) 
Simplifying and collecting like terms, we have 
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m; - Fe i-e) + (5-32) y tagen + (ZX) hm, 


2 10 2 10 90 90 45 
(i-r NLIS GF) (i-i) ) 

g. 35]. Hu eue g m 2 ee. eap prs 
similarly 


1 I , f 
m =f [x + shy + Az1T (V1) + Az2T (V2) + A23T (V3), y' + az21T' (v1) + a22T' (V2) 
+A23T'(v3)] 


1 5 15 OON 2|, x 
j=1 j=1 


5 v5\/,< 
+ TE h| y + ) az;hmj |, y' + azıhmı + a?;hm; + az3hm; 
j 


=1 


1 1 7 X415 1 7 v15 
m; = f (estre + (eS) wm, rarm ) rms, 


2 2 144° 72 36 144 72 
a ied eae 
Y Elgg" gp | rEg ma ag 24) "3 


1 v15 1 T 
mz =f |x+ 2 + ETE h,y + azıT (4) + az2T (v2) + az3T (v3), y' + azı T' (v4) + a32 
a3?T' (v5) + az3T'(v3)) 


1 VI5 l 3 l 3 
—fix- stip hyt+az3,h\y +). ayjhmy, + azzh| y +) ahm 
j=1 j=1 


3 


3 
+azzh | y' + 2. a4jhm; |, y" + » a5;hm; 
ja 


jzi 


Substituting for a;; from table 1 and simplify, we have 


m; =f (er (G42) a (oe) + (Z +E) nem, 


2 10 2° 10 9° 36 90 45 
ic EN p oum OR ) 
go ^ ^» X Eag 39] "4 Ng” 15 [Hz gg His 


The general solution of (2.02) is define as 
Yn+1 = Yn + baT (v1)  b;T(v;) + b3T (v3) 
5 r, 2, 415 5 15 4 ; 5 , /15 
Yn +ŽA(y + Am, + (24-33) hm; + (3; - X5) hms) +2n(y + (3 4 3) hm, + 


30 36 24 


“hm, J (s E) hms) +—h(y’ $ ($38) hm, + (2+) hm; tx hm) 


36 24 18 
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(2.05) 
Y'n+1 = Y'a + b4T'(v4) + b;T'(v;) + bzT' (v3) 
5 5 
= yu T ig T 1g" + 7g hs 
Simplify (2.05), the general solution for solving (2.02) is 
h? 
Yn41 = Yn + hy' + 36 [(5 + v15)m, + 8m, + (5 Lm v15)m;] 


5 
Yn = Y'a + —h[5(m, + m3) + 8m;] 


18 
where 
m,-fix-* ve h,y+ AT, y qe uan 3 NT h?m, 
2 10 : 2 10 90 90 45 
p IBN ed 2 wv15 5 V15 
+ 9 36 hmg, y +zghm + 9 15 hm; + 36 30 hm3 
= d ino (a 152 ES. 7 BY 9 ' 
m =f (x ey hy t (tv) m zh m * (4 7; ) ms y + 
($+ E) am hm; + (= - 2) hm) 
= 1 vi5 1 vi5 1 1 vi5 2 7 vi5 2 
m =f (rrt) +S +o) htm t (Gt Ge) heme + 


~h2m3,y' + (= + Hz) hm, + (E+ =) hm, + = hms) 


3.0 Analysis of the scheme 
The method for the second order ODEs is summarized in table 2: 
Table 2: Butcher’s tableau for second order ODEs 























C A A' 
1 v15 5 2 vi5 5 y/15 1 7 915 1 V15 
2 10 36 9 15 36 30 90 90 45 9 36 
1 5  V15 2 5 v15 7 , V15 1 7 V5 
2 36 24 9 36 24 144 72 36 144 72 
i v15 AA VI5\ (2 y V15 5 1 vi5 n V15 1 
2^ 10 36 30 9° 15 36 9^ 36 90 45 90 
1 8 5 1 8 1 
bT p 





=C | AXA  -R(T)e9e? 
bTXbT 
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Consistency: Table 1since the (Domain of T) is consistent, we concluded that table 2, is also 
consistent because it is the range of T. Since the Dim(T) is isomorphic to Rang(T) and T 
preserve the algebraic structure of Domain onto Rang(T ) 
Stability: The test stability equation is 
y'-Ay A«0 
For stability region, we set Z — Ah? 
The stability function is defined as R(z) where 
R(z) = I + (Zb? + Z?b'') (IZA — Z?A') te 
e = (,1,1)7, he (0,1), D = AxA is the coefficient matrix, I is the identity 3 x 
3 matrix, Z = Ah*,B = (bb, b3), b' = (b'4, b'5, b'4) weights. 
The A-stability Region is 
(R(z)) = (z: R(z) < 0 and |R(z)| < 1) 
Our method is automatically A-stable since the Dim(T) based on Gauss quadrature method 
which is A-stable (see Butcher [4]), and T preserve its algebraic structure onto its Range. 
Error: Since every Runge-kutta solution agree with Taylor’s series expansion of order p = 
6 
The error is C(h?*1) = C(h’) 
Proposition 1 (Error estimation formula) 


h 
If y» and yl) are approximate solutions of a Runge-kutta type method of order P with 


h 
step size h and 7 respectively, then the error C(h?*?) is 


h 
277 L3. m 
T 20*1.— 1 Yn+1 ^ Yn+1 
Proof 
h) (2) 
Let y,,,, and y,7, be approximation solution of Runge-kutta method of order p respectively. 


These solutions can be expanded into Taylor's series and Runge-kutta solutions agrees 
with the Taylor's series expansion up to the p terms and converges to exact solution 
Yn+1 -Thus we can write 


Yrs = Yep + CCPH) + Ry (x) ~ ou (1) (2.06) 
h pti 
Wai G) +C (5) + On (X) ... .. (Hi) (2.07) 


as p > oo, Rpa (x)& Qa (x) > 0, hence R (x) & Q,(x) can be ignored. Then subtracting 
(2.07) from (2.06) we have 


h pt+1 _ 
BORO _ (3) 27^ -1 
0— yu 7 Ynga + COP) — 
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2p+1 h x 
Cnt) = 2»*1—1 (3 = x) 


Since C(h?**) is the first neglect term of the series, our truncated error is C(h?*+). Thus 
our error is 


"T ien for rim 
C(hP**) = E, = 2P*1—4 Yn+1 — Yn+1 


Remark: Since our method is of order 6. The error of our Runge-kutta method is 


7 h 
E, = XL (na — x) (error estimation formula) 








(2.08) 


4.0 Numerical Experiments 
In this section we test the performances of our schemes by considering some problems with exact 
solutions to check quality, stability and implementation cost of our schemes 


Example 1 
1 
y'-xQo0,  y(00-1Ly()-; 
i i - ig foe aba 
Analytic solution: y(x) = 1+ ;In (=) 
Example 2 


y" —3y'+2y=x, yO)=1,y'(0) =0 
Analytic solution: y(x) = e* — se + zx + : 
Example 3 

y" = —100y, y(0) = 1,y'(0) 2 10 ,h = 0.01 
Analytic solution: y(x) = Cos (10 x) + Sin (10 x) 


Table 3: Comparison of the theoretical and approximate solutions for example 1 























x y(x) 3 à» ExactErrors | Approximation 
Ynti Yi) — Yn+i | errors 
(2.08) 
0.1 | 1.050041729278490 | 1.050041729277550 | 1.050041729278496 | 9.40 E (-13) 9.37 E (-13) 
0.2 | 1.100335347731070 | 1.100335347729130 | 1.100335347731050 | 1.94 E (-12) 1.94 E (-12) 
0.3| 1.151140435936470 | 1.151140435933350 | 1.151140435936440 | 3.12E (-12) 3.11 E (-12) 
0.4 | 1.202732554054080 | 1.202732554049490 | 1.202732554054060 | 4.59 E (-12) 4.60 E (-12) 
0.5 | 1.255412811883000 | 1.255412811876440 | 1.255412811882980 | 6.56 E (-12) 6.59 E (-12) 
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Table 4: Comparison of the theoretical and approximate solutions for example 2 

x y(x) 3 (5) ExactErrors | Approximation 

Ynti y(%i) — Yn+i | errors 
(2.08) 

0.1 | 0.989118849455522 | 0.989118849340112| 0.989118849453719| 1.15 E (-10) | 115 E (-10) 

0.2 | 0.952534234929220 | 0.952534234647043 | 0.952534234924814| 2.82E (-10) | 2.80 E (-10) 

0.3 | 0.883269707283120 | 0.883269706765726 | 0.883269707275048| 5.17 E(-10) | 5.13 E (-10) 

0.4 | 0.772669001271920 | 0.772669000428699 | 0.772669001258759| 8.43 E(-10) | 8.37 E (-10) 

0.5 | 0.610009899355840 | 0.610009898067594 | 0.610009899335757| 1.29 E(-09) | 1.28 E (-09) 

















Table 5: Comparison ofthe theoretical and approximate solutions for example 3 
































x y(x) yop Error [15] Exact Error 

0.01: 1.094837581924850 1.09483758192396 1.0 E (-08) 8.90E (-13) 
0.02: 1.178735908636300 1.17873590836347 2.3 E (-08) 1.55 E (-12) 
0.03} 1.25085669578695(0 1.25085669578490| 4.0E (-08) 1.98 E (-12) 
0.04. 1.31047933631154( | 1.31047933630941| 5.1 E(-08) | 2.13 E(-12) 
0.05} 1.35700810049458(0 1.35700810049258| 6.9 E(-08) | 2.00E (-12) 








5.0 Discussion of Results 





We observed from the three problems tested the approximate error method (2.08) 
converges to the exact errors. This shows that our method is good and can be used to solve 
accurately problems without exact solutions. (see tables 3,4 and 5). The new method has 
only three function evaluations m4,m;,m which are directly used for integration of the 
problems. The method is much better than that of [15] (see table 5) in terms of accuracy 
and implementation cost. 


6.0 Conclusion 

We have used three Gauss-quadrature points to derived a 6t^ order implicit Runge-kutta 
method for direct integration of second order differential equations. The method is A- 
stable, highly efficient and has simple coefficients with low implementation cost. Alongside 
we develop an accurate error formula for step size control. 
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